

In [79]:

    
%matplotlib inline
import numpy as np



In [75]:

    
class __Atom__:
    """ Crappy class of Atom for holding atomic
        symbol, charge and coordinates.
    """
    def __init__(self, Symbol, Coordinates):
        """ Label is atom symbol string,
            Coordinates is 3-tuple list
            
            Charges are taken from this list:
            https://en.wikipedia.org/wiki/Effective_nuclear_charge#Values
        """
        ChargeTable = {"C": 5.673,
                       "B": 4.680,
                       "H": 1.,
                       "O"
                      }
        self.Symbol = Symbol
        self.Coordinates = np.array(Coordinates)
        self.Charge = ChargeTable[Symbol]

    def GetCoordinates(self):
        return self.Coordinates
    
    def GetSymbol(self):
        return self.Symbol
    
    def GetCharge(self):
        return self.Charge



In [58]:

    
def CMDiagonal(Charge):
    """ Diagonal elements only depend on charge of one atom """
    return 0.5 * Charge**(2.4)

def CMOffDiagonal(Atom1, Atom2):
    """ Off-diagonal elements of Coulomb Matrix from Hansen et al. JCTC, 2013 """
    Denominator = np.sqrt(np.sum((Atom1.GetCoordinates() - Atom2.GetCoordinates())**2.))
    Numerator = Atom1.GetCharge() * Atom2.GetCharge()
    return Numerator / Denominator

def CalculateCoulombMatrix(Molecule):
    """ For a given dictionary input of atom classes,
        return the associated Coulomb matrix
    """
    NAtoms = len(Molecule)          # Number of atoms in the molecule
    CoulombMatrix = np.zeros((NAtoms, NAtoms), dtype=float)   # initialise array
    for i in xrange(NAtoms):        # Loop over atoms i
        for j in xrange(NAtoms):    # Loop over atoms j
            if i == j:              # Calculate diagonal element as effective charge
                CoulombMatrix[i, j] = CMDiagonal(Molecule[i].GetCharge())
            else:                   # Calculate off-diagonal elements
                CoulombMatrix[i, j] = CMOffDiagonal(Molecule[i], Molecule[j])
    return CoulombMatrix

Atoms

Atoms are defined as a list of strings here.



In [26]:

    
Atoms = ["C", "B", "H"]

Generate some random coordinates

Making some random $x,y,z$ coordinates up



In [27]:

    
Coordinates = []
for AtomNumber in range(len(Atoms)):
    Coordinates.append(np.random.rand(3))

Molecule

Here we've defined a molecule as a dictionary of atom objects.

We'll loop over the number of atoms, and initialise each dictionary key with an atom.

In an actual run, we'll need to parse information from an xyz file somehow instead.



In [76]:

    
Molecule = dict()



In [77]:

    
for Index, Atom in enumerate(Atoms):
    Molecule[Index] = __Atom__(Atom, Coordinates[Index])

Calculate the Coulomb matrix

Call our function, using the Molecule dictionary as input to calculate a Coulomb matrix.



In [81]:

    
CM = CalculateCoulombMatrix(Molecule)



In [ ]:

    
def ReadXYZ(File):
    f = open(File, "r")
    fc = f.readlines()
    f.close()
    NAtoms = len(fc)
    Coordinates = []
    for line in range(NAtoms):
        Coordinates.append(fc[line].split())
    return Coordinates